# Import necessary packages
from IPython.display import HTML, Image
import numpy as np
# import plotting packages
import matplotlib.pyplot as plt
import plotly.graph_objects as go
If at time $t=0$, $I_0 = 5$ and then $I_1 = 15$. We can calculate the rate of infection by the following
$$ I_1 = I_0 + x S_0 I_0 $$Solve for $x$ and use this as the new rate!
Up to now, we have been focusing on deterministic modeling. In case where we cannot explain a system analytically or data collected directly, we can simulate the behavior.
Processes with an element of chance involved are called probabilistic, as opposed to deterministic, processes.
Monte Carlo methods include a large class of computational algorithms that rely on repeated random sampling to obtain numerical results.
Thus, these methods are probabilistic.
Let us consider the equation of a circle
$$ x^2 + y^2 = 1 $$and we wanted to calculate the area of the circle (without knowing the formula).
def area_circle(x,y,n):
# initiate the count for inside the circle
in_c = 0
# check if the point (x,y) is inside the circle
for i in range(n):
if x[i]**2 + y[i]**2 <= 1:
# increment the counter by 1
# increase the number of points in the circle by 1
in_c = in_c + 1
#print('hey, (x,y) is in the circle')
# calculate the fraction of points inside the circle
in_cfrac = in_c / n
# calculate the area of the quarter circle
areac = in_cfrac
return areac
# determine number of points
n = 10_000_000
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=0, high=1, size=n)
y = np.random.uniform(low=0, high=1, size=n)
x.shape
(10000000,)
area_circle(x,y,n)*4
3.14175
# create/allocate memory for area
my_area = np.zeros(4)
# fill each element with the calculated area
my_area[0] = area_circle(x,y,5)
my_area[1] = area_circle(x,y,10)
my_area[2] = area_circle(x,y,100)
my_area[3] = area_circle(x,y,1000)
fig = go.Figure()
fig.add_trace(go.Scatter(x=[5,10,100,1000], y = my_area, mode='markers'))
We can do this for more than circles! If we wish to calculate the area under the curve of $y = f(x)$ over $ a \le x \le b$ by constructing a bounding rectangle with height $M$ and base length $b-a$.
Then,
$$ \frac{\text{area under curve}}{\text{area of rectangle}} \approx \frac{\text{number of points counted below curve}}{\text{total number of random points}} $$Input: Total number $n$ of random points to be generated in the simulation.
Output: Approximate area under $y=f(x)$ over $a \le x \le b$
1) Initialize: COUNTER = 0.
2) For $i = 1,2,\ldots,n$, do Steps 3–5.
3) Calculate random coordinates $x_i$ and $y_i$ that satisfy $a \le x_i \le b$ and $0\le y_i \le M$.
4) Calculate $f(x_i)$ for the random $x_i$ coordinate.
5) If $y_i \le f(x_i)$, then increment the COUNTER by 1. Otherwise, leave COUNTER as is.
6) Calculate AREA = $M(b - a)$ COUNTER/ n.
7) OUTPUT (AREA)
Use the Monte Carlo Area Algorithm to approximate the area under the curve of $y = \cos(x)$ from $-\pi/2 \le x \le \pi /2$.
np.pi
3.141592653589793
# np.cos() -- cosine
# np.pi -- pi
fig = go.Figure()
fig.add_trace(go.Scatter(x = x, y = np.cos(x), mode='markers'))
fig.show()
def area_cos(x,y,n):
# initiate the count for inside the circle
in_c = 0
# check if the point (x,y) is below f(x)
for i in range(n):
if y[i] < np.cos(x[i]):
# increment the counter by 1
in_c = in_c + 1
# calculate the fraction of points below f(x)
in_cfrac = in_c / n
# calculate the area
areac = 1 * np.pi * in_cfrac #[M]*(b-a)*in_cfrac
return areac
# determine number of points
n = 10
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=-np.pi/2, high=np.pi/2, size=n)
y = np.random.uniform(low=0, high=1, size=n)
area_cos(x,y,n)
2.5132741228718345
Use the Monte Carlo Area Algorithm to approximate the area under the curve of $y = x^3 - \sin(x)$.
def area_sin(x,y,n):
# initiate the count for inside the circle
in_c = 0
# check if the point (x,y) is below f(x)
for i in range(n):
if y[i] < (x[i]**3) - np.sin(x[i]):
# increment the counter by 1
in_c = in_c + 1
# calculate the fraction of points below f(x)
in_cfrac = in_c / n
# calculate the area
areac_new = 1 * np.pi * in_cfrac #[M]*(b-a)*in_cfrac
return areac_new
n = 10
# simulate x,y points (each point is inside [-1,1])
x = np.random.uniform(low=-0.5, high=0, size=n)
y = np.random.uniform(low=0, high=1, size=n)
area_sin(x,y,n)
0.9424777960769379